clear all

x = [0.1;0;0];
n = numel(x);
tau = 5;
v  = .5;
f = inline(sprintf('[x(2)*%0.2f;-x(1)*%0.2f;(x(2)-x(3))/%0.2f]',v,v,tau),'x');


N = 1000;

X = zeros(n,N+1);
X(:,1) = x;
dt = 0.05;
%J = [0 1 ; -1 0];
%J = [0 -1 ; 1 0]';
J = [0 -1*v 0; 1*v 0 1/tau; 0 0 -1/tau]';

r = 0* rand(n,N);

for i = 2:N+1
     X(:,i) =  X(:,i-1)+linsolve(J,(expm(J*dt)-eye(n))*f(X(:,i-1))) + r(:,i-1);
end


figure
plot(X')